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! Abstract 

We study the bound states of two-dimensional bright solitons in nonlocal nonlin- 
ear media. The general properties and stability of these multisolitary structures are 
investigated analytically and numerically. We have found that a steady bound state 
p' ■ of coherent nonrotating and rotating solitary structures (azimuthons) can exist above 



some threshold power. A dipolar nonrotating multisoliton occurs to be stable within 
the finite range of the beam power. Azimuthons turn out to be stable if the beam 
power exceeds some threshold value. The bound states of three or four nonrotating 
solitons appear to be unstable. 



The recent experimental observations of spatial solitons in nonlocal media such as nematic 
liquid crystals [1,2], lead glasses [3], renewed an interest to coherent structures in spatially 
nonlocal nonlinear media. In the spatially nonlocal media the nonlinear response depends 
on the wave packet intensity at some extensive spatial domain. Nonlocality is a key feature 
of many nonlinear media. It naturally appears in different physical systems such as plasmas 
[4-6], Bose- Einstein condensates [7], optical media [8], liquid crystals [1], and soft matter 
[9]. Different types of single (solitons , vortices) and composite (soliton clusters) solitary 
structures have been predicted and experimentally observed in various nonlinear media [10]. 
Rotating scalar multipoles (azimuthons) were first introduced in Ref. [11] for models with 
' local nonlinearity. 

A composite soliton structure, or a multi-soliton complex is a self-localized state which is 
a nonlinear superposition of several fundamental solitons. Stationary multi-soliton structures 
in nonlocal media were considered first in Refs. [12,13]. However, a question of the stability 
of the multisolitons is still an open problem. In particular, the form of the nonlinear nonlocal 
response is crucial for the multisoliton stability. Stable two-dimensional (2D) rotating [14] 
and nonrotating [15] dipole solitons in a medium with Gaussian response function have been 
numerically observed. Authors of Refs. [16] investigated a more realistic model with the 
Helmholtz type response function (often referred to as the model with thermal nonlinearity) 
and showed that multipole vector solitons can be made stable. This model describes, in 
particular, the nonlinear nonlocal response of some thermo-optical materials, liquid crystals, 
and partially ionized plasmas [1,3,5]. Very recently experimental observation of scalar 
multipole solitons in the optical media with thermal nonlocal nonlinearity has been reported 
[17]. It is particularly remarkable that all multisolitons investigated in Ref. [17] are found to 
be unstable, but with a large parameters range where the instability is weak, so that their 
experimental observation is possible. In this Letter, using direct 2D simulations for a scalar 
model with thermal nonlocal nonlinearity, we find a class of radially asymmetric nonrotating 
two-dimensional soliton solutions and show that, at some input power, the dipole-mode 



scalar solutions are stable, while the tripoles and quadrupoles are always unstable, but can 
survive over quite considerable distances. The stability window for dipolar multisolitons is 
found. We find also rotating dipole and quadrupole solutions (azimuthons) which turn out 
to be stable if the beam power exceeds some threshold value. 

The basic system of equations, written in appropriate dimensionless variables, is 

i-z- + A_l# + 0* = O, (1) 
az 

a 2 9 - A ± 9 = |^| 2 , (2) 

where Aj_ = d 2 / dx 2 + d 2 / dy 2 is the transverse Laplacian. Equations (0) and (J2J) describe the 
propagation in z-direction of the electric field envelope *$?(x, y, z) coupled to the temperature 
perturbation 8(x,y,z) in a plasma [5]. The identical model describes the wave field \P and 
spatial distribution of the molecular director 9 in nematic liquid crystals [1]. 

The parameter a stands for the degree of nonlocality of the nonlinear response. In the 
limit a 2 3> 1, Eqs. ((H) and (J2J reduce to the usual nonlinear Schrodinger equation; the 
opposite case « 2 C 1 corresponds to a strongly nonlocal regime. For plasmas parameter 
a characterizes the relative energy that electron with mass m delivers to a heavy particle 
with mass M during single collision (a 2 ~ 2m/ M). Therefore, the thermal self-focusing in 
a plasma occurs in a strongly nonlocal regime. The model used in Ref. [3] to describe the 
propagation of the wave beam in optical media with thermal nonlinearity corresponds to the 
specific nonlocal limit a — > 0. In nematic liquid crystals, the degree of nonlocality can be 
modulated by changing the pretilt angle 9 of molecules through bias voltage V [18]. As V 
increases, the degree of nonlocality decreases, so that the degree of nonlocality can be tuned 
from local regime to a strongly nonlocal one. 

Equations (JIJ and (J2J) can be rewritten as a single integro-differential equation 



dz 



+ A ± * + ^ J R(\r-r'\)\^(r')\ 2 dr' = 0, (3) 



with the kernel R(£) = Ko(a£) / (2tt) , where K (z) is the modified Bessel function of the 
second kind of order zero. It is seen, that the nonlinearity in Eq. (j3J) has essentially nonlocal 
character. 

We look for stationary solutions of Eqs. (0) and (J2J) in the form ^(x, y, z) = if)(x, y) exp(i\z) 
so that i/;(x, y) obeys the equations 

- \ip + A ± ip + 9ip = 0, (4) 

a 2 9- A ± 9 = |^| 2 , (5) 

and we do not assume the radial symmetry of ip(x,y). For the numerical modeling we use 
the rescaled variables ip — ► ip/a, 9 — > 9/ a 2 , z — > za 2 , (x,y) — > (x, y)a, X — > A/a 2 , so 
that a = 1 in Eqs. (j3J and ©. Note, a strongly nonlocal regime (a 1) corresponds to 
large values of the scaled propagation constant A. Imposing periodic boundary conditions 
on Cartesian grid and choosing an appropriate initial guess we have found a class of radially 
asymmetric multipole localized solutions by using the relaxation technique similar to one 
described in Ref. [19]. The real (or containing only a constant complex factor) function 
ip(x,y) corresponds to nonrotating solitary structures. Examples of such nonrotating mul- 
tipole solitons, namely, a dipole, a tripole, and a quadrupole are presented in Fig. [T] The 
nonrotating multipoles consist of several fundamental solitons (monopoles) with opposite 



phases. The complex function ip(x,y) with a spatially modulated phase corresponds to ro- 
tating structures (azimuthons) with nonzero angular momentum. In Fig. |2] we demonstrate 
two examples of the azimuthons for the nonlocal model described by Eqs. (J!Q) and 

To gain the better insight into the properties of the stationary nonrotating multisolitons 
we have performed the analytical variational analysis. A stationary nonrotating multisoliton 
can be described by the trial function of the form: 

*{x,y,z) = hf(x,y)e i >«, (6) 

where the shape of the dipolar multisoliton can be approximated by the antisymmetric 
superposition of two gaussian functions: 

ffa v ) = e-sM^- rf / 2 ) 2 } - e~^i y2+{x+d/2)2 }, 

where a and d characterize the radius of each monopoles and the distance between them, 
respectively. For the bound state of S identical single solitons with the centers at (x s ,y s ), 
where s = 1 . . . S, we have used the trial function of the form: 

s=l 

where (£,77) = (x,y)/a, g s = ±1 gives the phase of the s-th soliton. 

As known, the steady-state corresponds to the stationary point of the Hamiltonian 

H = f (|V^| 2 - l-9\y\ 2 }d 2 r (7) 



2 



N= / \V\ 2 d 2 r. (8) 



at the fixed number of quanta iV 



In the variational approach, Hamiltonian if (a, b) is the function of two variational parameters 
a and b = d/a, where a characterizes the width of each monopole and d is the distance 
between two neighboring out-of-phase solitons. The third parameter of the trial function 
(jHJ), the amplitude h, has been eliminated in terms of a, b, and N using the normalization 
condition (JHJ). 

Results of the variational analysis are found to be in very good agreement with our 
numerical simulations for stationary nonrotating multisolitons. Figure El (a) shows the beam 
power N versus the propagation constant A for monopoles, dipoles, tripoles, and quadrupoles. 
The threshold power for existence of the stationary bound state of S solitons is very close 
to the threshold for formation of S single unbound fundamental solitons. Moreover, the 
"binding energy" Ns — SNi (where N% is the power for single fundamental soliton) is very 
small for A < 1.9, thus these multisolitons should readily decay into inbound solitons. As 
is seen from Fig. |3] (c), the distance between neighboring out-of-phase solitons d increases 
dramatically when A < 1.9 which lead to decreasing of interaction between monopoles. 

We next addressed the stability of these multipole solutions and study the evolution 
(propagation) of the dipoles, tripoles, quadrupoles, and azimuthons in the presence of small 
initial perturbations. We have undertaken extensive numerical modeling of Eqs. (JIJ and (J2J) 
initialized with our numerically computed multipole solutions with added gaussian noise. 
The initial condition was taken in the form ip(x, y) [1 +£$(x, y)}, where ip(x, y) is the numer- 
ically calculated exact multipole solution, <&(x,y) is the white gaussian noise with variance 
a 2 = 1 and the parameter of perturbation e = 0.005-^0. 1 . Spatial discretization was based on 
the pseudospectral method and "temporal" z-discretization included the split-step scheme. 



Depending on the parameter A, we observed three different regimes of the nonrotating 
dipole propagation, which are presented in Fig. 0] (for e = 0.02). The first regime corresponds 
to the region A < A cr , and we found A cr ~ 1.9. If A < A cr , the initial dipole splits in two 
monopoles which move in the opposite directions without changing their shape and without 
radiation, i. e. the monopoles just go away at infinity. This type of the evolution is shown in 
Fig-Efa). The splitting of the dipole into two monopoles can be easily understood. It is seen 
from Fig. HJthat the bound energy 5N = Ndi P — 2N mon in the dipole tends to almost zero as 
A approaches A cr ~ 1.9. This explains why the dipole with A < A cr can be easily (i. e. under 
the action of extremely small initial perturbations) split into two monopole-type solitons. A 
similar behavior, i.e. the decay of the initial dipole into two stable moving monopoles below 
some critical value of A, was observed for the model with a Gaussian response function R(C.) 
in Eq. © [15]. 

The second regime corresponds to the region A cr < A < Xth, where X t h ~ 4. The numerical 
simulations clearly show that in this range of the parameter A the dipoles are stable with 
respect to initial noisy perturbations. If the parameter of perturbation e is not too large, 
the dipoles survive over huge distances (z > 3000). The stable propagation of the dipole is 
illustrated in Figs. Efb) (for A = 3.5 and e = 0.02). 

The further (after A^ ~ 4) increasing the parameter A sharply shortens the propagation 
distances at which the dipole survives, and, the dipoles with A > A^ are unstable. The 
typical decay of the unstable dipole above the threshold value A^ of the rescaled propagation 
constant is shown in Figs. Efc). Thus, the stable dipoles exist only within a finite, rather 
narrow range of the propagation constants A. 

Figure El illustrates the propagation of the tripole and quadrupole for A = 2, i.e. in the 
region, where the dipole is stable. Generally speaking, the tripoles and quadrupoles turn out 
to be unstable, but for A cr < A < \ t h they can survive on the quite considerable (compared 
to the characteristic diffraction length) distances and, thus, can be experimentally observed. 
Tripoles and quadrupoles with A < A cr decay into three and four monopoles respectively. 

We observed two different scenarios for the azimuthon propagation. If A (i. e. the beam 
power) is small enough (A < 15), the azimuthons are unstable and split into fundamental 
solitons (monopoles) which move away from each other. An example of such splitting for the 
azimuthon with two intensity peaks and A = 10 is presented in Fig.|^a). For sufficiently large 
A (A > 15), the azimuthons turn out to be stable. Stable propagation of the azimuthons 
with two (rotating dipole) and four (rotating quadrupole) intensity peaks and A = 80 is 
shown in Figs. EKb),(c). Note that stable rotating dipoles were also observed in Ref. [20]. 
Numerically estimated rotational velocity of the stable azimuthon in Fig. |^b) is uo = 1.15 
so that it survives over hundreds rotational periods. 

The dynamics of the nonrotating dipole in our model Eqs. and (J2J) is in sharp 
contrast to the dipole propagation in the model with a Gaussian response function -R(£) in 
Eq. (J2I), where the stable nonrotating dipoles were observed for all A > A^, where \ t h is 
some threshold value [14,15]. A qualitatively different behavior of the dipoles in these models 
seems to be related to the regularity properties of the functions R(£) in Eq. (jSJ) - for the 
model described by Eqs. (JTJ) and (0) the function R(£) has a singularity at zero. The strong 
dependence of a stability criteria for multisoliton solutions on the regularity properties of 
the kernel R in Eq. (J3J) was discussed also in Ref. [15,20]. The corresponding theoretical 
problem seems to be intriguing and highly nontrivial. 

In conclusion, we have studied the bound states of the out-of-phase nonrotating two- 
dimensional solitons and rotating multisolitons in nonlocal nonlinear media. We have demon- 
strated that stationary nonrotating dipolar, tripolar, quadrupolar multisolitons and rotating 
multisolitons (azimuthons) may exist in media with thermal nonlocal nonlinearity, if the 



beam power is above some critical value. We have investigated stability of the multisolitons 
using direct numerical simulations. The nonrotating tripoles and quadrupoles are found to 
be unstable with respect to decay into fundamental solitons or fusion into single soliton, de- 
pending on the beam power. At the same time, we have observed robust propagation of the 
azimuthons with sufficiently high energy and a dipolar multisoliton with moderate energy. 
Thus, our theoretical predictions open the prospects for the experimental observations of 
the stable azimuthons and the bound states of two out-of-phase wave beams in media with 
thermal self-focusing nonlocal nonlinearity. 
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Figure 1: Numerically found stationary localized nonrotating solutions of Eqs. (jU and (J5J): 
(a) monopole with A = 10; (b) dipole with A = 10; (c) tripole with A = 15; (d) quadrupole 
with A = 20. 




Figure 2: Two examples of the azimuthons with A = 4.5. Intensity |^| is shown in the (x, y) 
plane, a) solution with two intensity peaks; b) solution with four intensity peaks. 




Figure 3: (a) Number of quanta iV vs propagation constant A. Solid curves for the bound 
states of S out-of-phase solitons, dashed curves for SNi (variational results). The integers 
near the curves indicate S, the number of solitons. 
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Figure 4: (a) Splitting of the dipole with A = 1.5 into two monopoles; (b) stable propagation 
of the dipole with A = 3.5; (c) unstable dipole with A = 10. 




Figure 6: (a) Splitting of the azimuthon with two intensity peaks and A = 10; (b) stable 
propagation of the azimuthon with two intensity peaks and A = 80; (c) stable propagation 
of the azimuthon with four intensity peaks and A = 80. 



